Python implementation of Kaiser windowed sinc functions

In [1]:

    

%pylab inline
from scipy.special import i0 as bessi0


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

x = linspace(-4, 4, 51)

taper = sqrt(1 - (x/4)**2)**0.5
kaiser = i0(6.31 * taper) / i0(6.31)
sincx = sinc(x)
kws = sincx * kaiser

plot(x, taper, label='Original window')
plot(x, kaiser, label='Kaiser window')
plot(x, sincx, label='Normalized sinc function')
plot(x, kws, label='Kaiser windowed sinc function')
legend(loc='upper left', bbox_to_anchor=(1,1))


    

    
Out[2]:





<matplotlib.legend.Legend at 0x10d80bf50>






    

In [3]:

    

Zi, Xi = np.mgrid[0:9,0:9]


    

In [4]:

    

Xi


    

    
Out[4]:





array([[0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8]])

In [5]:

    

HC_KAISER = {
    1:  1.24,
    2:  2.94,
    3:  4.53,
    4:  6.31,
    5:  7.91,
    6:  9.42,
    7:  10.95,
    8:  12.53,
    9:  14.09,
    10: 14.18,
}

def KaiserWindowedSinc(ireg, offset):
    '''
    Finds 2D source terms to approximate a band-limited point source, based on

    Hicks, Graham J. (2002) Arbitrary source and receiver positioning in finite-difference
        schemes using Kaiser windowed sinc functions. Geophysics (67) 1, 156-166.

    KaiserWindowedSince(ireg, offset) --> 2D ndarray of size (2*ireg+1, 2*ireg+1)
    Input offset is the 2D offsets in fractional gridpoints between the source location and
    the nearest node on the modelling grid.
    '''
     
    from scipy.special import i0 as bessi0
    import warnings

    ireg = int(ireg)
    try:
        b = HC_KAISER.get(ireg)
    except KeyError:
        print('Kaiser windowed sinc function not implemented for half-width of %d!'%(ireg,))
        raise

    freg = 2*ireg+1

    xOffset, zOffset = offset

    # Grid from 0 to freg-1
    Zi, Xi = numpy.mgrid[:freg,:freg] 

    # Distances from source point
    dZi = (zOffset + ireg - Zi)
    dXi = (xOffset + ireg - Xi)

    # Taper terms for decay function
    with warnings.catch_warnings():
        warnings.simplefilter('ignore')
        tZi = numpy.nan_to_num(numpy.sqrt(1 - (dZi / ireg)**2))
        tXi = numpy.nan_to_num(numpy.sqrt(1 - (dXi / ireg)**2))
        tZi[tZi == inf] = 0
        tXi[tXi == inf] = 0

    # Actual tapers for Kaiser window
    taperZ = bessi0(b*tZi) / bessi0(b)
    taperX = bessi0(b*tXi) / bessi0(b)

    # Windowed sinc responses in Z and X
    responseZ = numpy.sinc(dZi) * taperZ
    responseX = numpy.sinc(dXi) * taperX

    # Combined 2D source response
    result = responseX * responseZ

    return result


    

In [6]:

    

a = KaiserWindowedSinc(9, [0.5, 0.5])


    

In [7]:

    

imshow(a, interpolation='None', cmap=cm.gray_r)


    

    
Out[7]:





<matplotlib.image.AxesImage at 0x10dc718d0>






    

In [8]:

    

plot(a[9,:])


    

    
Out[8]:





[<matplotlib.lines.Line2D at 0x10de78c10>]






    

In [8]: